/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
    \file poly.hpp
    \brief An efficient implementation of a polynomial.

    This class include usefull and fast methods to evaluate polynomial,
    find polynomial roots, differentiate, integrate, data fitting and
    interpolation.
*/

# ifndef POLY_HPP
# define POLY_HPP 1

// A number which is less than POLYEPS is zero for us.
#ifndef POLYEPS
#define POLYEPS 1e-16
#endif

#include <iostream>

#include "typedefs.hpp"

namespace spectral {

/**
    \class Poly
    \brief Polinomial rappresentation class.

    Polynomials are written in the form:
    \f[
	\sum_{i=0}^n a_i x^i
    \f]
    so the only parameter that characterizes each polynomial
    is the coefficients' array \b a.
*/

/**
    \fn Poly::Poly()
    \brief Default polynomial constructor.

    This constructor will create a void polynomial of degree ``-1''
*/

/**
    \fn Poly::Poly(int const n)
    \brief Build an empty polynomial

    This constructor will initializate a zero polynomial with degree n.
    WARNING: as this will build an empty polynomial with a degree which is nonzero,
    working with this polynomial if the coefficients are not re-assigned could
    lead to strange behaviour
*/

/**
    \fn Poly::Poly(double const & x)
    \brief Convert a double to a polynomial.

    This constructor is used to convert a number into a polynomial. This
    is useful if we want to sum or multiply a double to a polynomial.
*/

/**
    \fn Poly::Poly(Poly const & p)
    \brief Poly copy constructor.

    This is the copy constructor for Poly
*/

/** \fn Poly & Poly::operator=(Poly const & p)
    \brief Copy assignment.
*/

/** \fn Poly::Poly(double const * const & coefficients, int degree)
    \brief Build a polynomial with given coefficients.

    This constructor will create a polynomial of degree \e degree :
    \f$ p(x) = a_i x^i \f$ where \f$ i \f$ goes from 0 to \e degree and
    \f$ a_i = \f$ \e coefficients[i] .
*/

/** \fn int Poly::degree() const
    \brief Returns the degree of the polynomial.

    The degree will be 0 for the constants and -1 for the void polynomial.

    This is an \b inline function.
*/

/** \fn double Poly::a(int n) const
    \brief Returns \f$ a_n \f$.

    This is an \b inline function.
*/

/** \fn double & Poly::operator[](int n)
    \brief Returns \f$ a_n \f$

    Please notice that this will return a <B>reference to</B> \f$ a_n \f$
    so the instruction \e p[3]=2; will be valid. For this reason the
    use of the const function .a is prefferred where possible.
*/

/**
    \fn void Poly::setdegree(int n)
    \brief Change the degree of the polynomial

    Warning. This will nuke the polynomial coefficients!

    This will change the degree of the polynomial to \e n. A new
    vector of coefficients is reserved in memory. The old one is
    deleted!
*/

/**
    \fn void Poly::setdegree()
    \brief Computes the new degree of the polynomial

    This will compute the new degree of the polynomial if the
    coefficients of the highest degree monomials are zero.
*/

/** \fn inline double Poly::operator()(double x) const
    \brief Evaluate the polynomial in \e x

    This function uses the Horner method to evaluate the polynomial.
    This method has linear complexity and doesn't evaluate any power of \e x
    so it is \b very fast!
*/

/** \fn Poly Poly::differentiate() const
    \brief Calculate the derivative of the polynomial

    This function outputs a new polinomial \f$ b(x) \f$
    such that: \f$ b(x)=\frac{d}{dx}p(x) \f$
*/

/** \fn Poly Poly::primitive() const
    \brief Calculate the primitive of the polynomial

    This function outputs a new polinomial \f$ b(x) \f$ such
    that: \f$ \frac{d}{dx} b(x) = p(x) \f$.
    The chosen arbitrary constant in \f$ b(x) \f$ is 0.
*/

/** \fn double Poly::integrate(double const & a, double const & b) const
    \brief Calculate the integral of the polynomial on a given intervall (a,b)

    This method does not call the primitive and evaluate method and will not
    instanziate a new polynomial object so if should be preferred over
    a joint use of the evaluate and primitive methods for efficiency reasons.
*/

/** \fn Poly & Poly::operator+=(Poly const &p)
    \brief Sum and assignment operator

    Please notice that if \f$ a(x) \in P_n \f$ and
    \f$ b(x) \in P_m \f$ doing \e a+=b or
    \e b+=a will not have the same complexity since if
    \f$ n > m \f$ \e b+=a will require
    to create a temporary vector of coefficients
*/

/** \fn Poly & Poly::operator*=(double const& alpha)
    \brief Multiplication by a scalar
*/

/** \fn Poly & Poly::operator-=(Poly const & p)
    \brief Difference and assignment operator

    This is designed to reuse the code from the sum and assignment operator
    without the need to build a new polynomial \e -*this
*/

/**
    \fn bool Poly::operator==(Poly const & p) const
    \brief true if the two polynomials are equal
*/

/**
   \fn Poly deflate(Poly const & p,double const & r)
   \brief Horner deflaction

   Given a root r of a polynomial p this function returns the deflacted
   polynomial q
 */

/**
 * \fn Poly operator-(Poly p)
 * \brief Unary Minus
 *
 * Returns -p
 */

/**
 * \fn Poly operator+(Poly a,Poly b)
 * \brief Sum Operator
 */

/**
 * \fn Poly operator-(Poly a, Poly b)
 * \brief Returns a-b
 */

/**
 * \fn Poly operator*(Poly a, double x)
 * \brief Returns a*x
 *
 * Overloading of the * operator for the multiplication of a Poly with a double.
 */

/**
 * \fn Poly operator*(double x, Poly a)
 * \brief Return x*a
 *
 * Overloading of the * operator for the multiplication of a double with a Poly.
 */

/**
 * \fn double newton(Poly const & p, double delta, double eps, double x0, int nmax)
 * \brief Classical Newton Method
 *
 * Computes the root of a given polynomial near the point x0
 * with a classical fixed-point method.
 */

/**
    \fn double * horner(Poly const & p, double delta, double eps, int nmax)
    \brief Newton-Horner Method

    Computes all the roots of a given polynomial using the Newton-Horner Method.
    The roots are supposed to be all __real__!
*/

/**
    \fn void legendre(Poly & L)
    \brief Legendre polynomial generator

    Computes the Legendre polynomial of degree \e degree over the
    reference interval [-1,1]. The given polynomial should be initializated with
    the correct degree:
    \code
    int n=10;
    Poly L(n);
    legendre(L);
    \endcode
    Will make L a Legendre polynomial of degree 10.
*/

/**
    \fn std::ostream & operator<<(std::ostream & strm,Poly P)
    \brief Overloaded << operator

    Used to print Polynomials
*/

/**
    \fn std::istream & operator>>(std::istream & strm,Poly & P)
    \brief Overloaded >> operator

    This is used to fill the polynomial coefficients with the input stream of numbers
    stranting from the lower orders monomials. Excessive input is thrown away withtout
    any check.
*/

class Poly {
    public :
        // Constructors and destructors:
        Poly();
        Poly(int const n);
        Poly(double const & x);
        Poly(Poly const & p);
        Poly(double const * const & coefficients, int degree);
        Poly& operator=(Poly const & p);
        ~Poly();
        // Interface:
        int degree() const;
        double a (int n) const;
        double & operator[](int n) const;
        void setdegree(int n);
        void setdegree();
        // Special operations:
        inline double operator()(double x) const;
        Poly differentiate() const;
        Poly primitive() const;
        double integrate(double const & a, double const & b) const;
        // Standard polynomial operations:
        Poly & operator+=(Poly const &p);
        Poly & operator*=(double const &alpha);
        Poly & operator-=(Poly const &p);
        bool operator==(Poly const & p) const;
        friend std::ostream & operator<<(std::ostream& strm,Poly P);
        friend std::istream & operator>>(std::istream& strm,Poly & P);
    private :
        int _P_degree;     /**< \brief Degree of the polynomial degree=-1 for void polynomial. */
        double * _P_a;     /**< \brief Polynomial coefficients */
};

Poly deflate(Poly const & p, double const & r);

Poly operator-(Poly p);

Poly operator+(Poly a, Poly b);
Poly operator-(Poly a, Poly b);

Poly operator*(Poly a, double x);
Poly operator*(double x, Poly a);

double newton(Poly const & p,
              double delta,
              double eps,
              double xnew,
              int nmax);   /**< Newton method for the polynomial */

double * horner(Poly const & p,
                double delta,
                double eps,
                int nmax=1000);

void legendre(Poly & L);

std::ostream & operator<<(std::ostream & strm,Poly P);

std::istream & operator>>(std::istream & strm,Poly & P);

// =======================================================
// End of Proper Header
// =======================================================

inline double Poly::operator()(double x) const {
    double u(_P_a[_P_degree]);
    for(int i(_P_degree-1);i>=0;--i){
        u=u*x+_P_a[i];
    }
    return u;
}

} // namespace
#endif
